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Using a LDG method for solving an 
inverse source problem of the 
time-fractional diffusion equation 


S. Yeganeh, R. Mokhtari* and S. Fouladi 


Abstract 


In this paper, we apply a local discontinuous Galerkin (LDG) method 
to solve some fractional inverse problems. In fact, we determine a time- 
dependent source term in an inverse problem of the time-fractional diffusion 
equation. The method is based on a finite difference scheme in time and a 
LDG method in space. A numerical stability theorem as well as an error 
estimate is provided. Finally, some numerical examples are tested to confirm 
theoretical results and to illustrate effectiveness of the method. It must be 
pointed out that proposed method generates stable and accurate numerical 
approximations without using any regularization methods which are neces- 
sary for other numerical methods for solving such ill-posed inverse problems. 


Keywords: Local discontinuous Galerkin method; Inverse source problem; 
Time-fractional diffusion equation. 


1 Introduction 


Recently, studying problems involving fractional order partial differential 
equations (PDEs) has attracted a lot of interests of scientists and engi- 
neers, see e.g. [1, 5-12, 15, 18-20, 22, 24-30, 33-35]. One of the fractional 
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PDEs which has been very popular is the fractional diffusion equation (FDE), 
see e.g. [1, 6-12, 15, 18-20, 22, 24-30, 33-35]. Nowadays instead of the classi- 
cal diffusion equations, FDEs have attracted wide attentions since a FDE 
is a generalization of a diffusion equation which can be used to describe an 
anomalous diffusion phenomenon such as super-diffusion or sub-diffusion. By 
replacing the standard time derivative with a time fractional derivative in a 
diffusion equation, a time-fractional diffusion equation is obtained. FDEs 
can be applied in modelling of some problems in porous flows, rheology and 
mechanical systems, models of a variety of biological processes, control and 
robotics, transport in fusion plasmas, and many other areas of applications. 
Analytical or numerical studying of the direct problems corresponding to 
the time-fractional diffusion equations have been carried out extensively in 
the recent years, see e.g. [6, 7,9, 11, 12, 15, 30] as well as references cited 
therein [27,28] for some subjects related to the obtaining some uniqueness 
and existence results, establishing maximum principle, finding analytical so- 
lutions, applying some numerical methods such as finite element methods or 
finite difference methods. 


If in an initial or initial-boundary value problem some parts of data such 
as boundary data, or initial data, or source term or even some coefficients of 
the main equation may not be given, we have encountered an inverse problem. 
Inverse problems are appeared in many practical situations, and for solving 
them we need to some additional measured data, see e.g. [1,8,10,16—25,27-29, 
32-35] and references cited therein. In some inverse problems, we need to find 
space-dependent source term [27,32] or time-dependent source term [28] or 
even space- and time-dependent source term [23]. In this paper, we focus our 
attention on finding the time-dependent source term in a fractional inverse 
problem which is one of the interesting and novel inverse problems. One of 
the pioneering scientists in this subject is Murio [18-20]. In the following 
we mention some of works which have been published after that. An inverse 
problem for determining the order of fractional derivative and the diffusion 
coefficient in a FDE has been considered in [1] where a uniqueness result has 
been also obtained. A backward problem for the time-fractional diffusion 
equation has been solved by a quasi-reversibility regularization method in 
[13]. In [34,35] some Cauchy problems based on the time-fractional diffusion 
equation have been investigated on a bounded domain and on a strip domain, 
respectively. Qian [22] has investigated a modified kernel method for solving 
an inverse fractional diffusion equation. An inverse source problem for a 
FDE has been solved in [33]. Rundell et al. [10] have been investigated some 
nonlinear fractional inverse problems. Recently, some interesting works have 
been carried out by Wei et al., see e.g. [27-29]. In [10] known theoretical 
results and computational techniques for FDEs have been provided. 


In [28] the following initial-boundary value problem for the time-fractional 
diffusion equation has been considered 
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Deu = Ure + f(x)p(t), O0< a <1,0<t<T, 


u(0, t) = ko(t), O<t<T, 4) 
u(1,t) = ky(t), O<t<T, 
u(x,0) = d(x), 0<a¢<l, 


where D? is the Caputo fractional derivative of order a, i.e., 


1 * Au(a,s) ds 
Oy = 1 2 
eu a. Os “ase ee 2) 


in which I is the Gamma function. For a = 1, we have D?'u = uz. We assume 
that ko, ky and ¢ are given functions. Problem (1) will be a direct (forward) 
problem whenever f and p are known functions. Based on the availability of 
f or p, there are some inverse problems. If p is known and f is unknown we 
need to an extra condition such as 


u(x, T) = g(x), O0<a<l. 


We have investigated numerical solution of this inverse problem in [82]. 
In another inverse problem, f is known but p is unknown and an over- 
determination condition such as 


u(ax*,t) = g(t), 0<t<T, (3) 


will be needed where x* € (0,1) is an interior measurement location. It must 
be pointed out that although these two inverse problems are very similar 
they are completely different. The inverse source problem which we consider 
here is to determine (p, u) based on problem (1) and condition (3). Using the 
main equation in (1) and Eq. (3), we have 
1 2u(ar,t 
v(t) = a5 (Dealt) - SP?) (4) 


We need (4) just for theoretical purposes. 


The above mentioned inverse source problem is an ill-posed problem [28]. 
Sakamoto et al. [24] have given a stability estimate for the inverse source 
problem and mentioned that the uniqueness of p is guaranteed if f 4 0, but 
they did not present any numerical method. This inverse source problem has 
been solved numerically by Wei et al. [28] using a regularized method based 
on the boundary element discretization for recovering a stable approximation 
to p. In [28], a numerical method has been applied to various examples in 
particular some examples with none-smooth data which lead to the none- 
smooth solutions. Therefore, we decided to apply the discontinuous Galerkin 
method to the above mentioned inverse source problem. To the best of our 
knowledge, for the inverse source problem with or without fractional opera- 
tors, the development of the DG methods remains limited and there are a few 
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results, see e.g. [32] where a space-dependent source term is determined in 
a time-fractional diffusion equation by using a local discontinuous Galerkin 
method. Of course, some DG methods have been applied successfully for the 
forward fractional diffusion equation. For example, Hesthaven et al. [6,30] 
have been solved some space-fractional diffusion equations using a local dis- 
continuous Galerkin method in a semi-discrete regime and Yang Liu et al. [14] 
have developed a LDG methos combined with WSGD approximation for a 
time fractional subdiffusion equation. 

In the present paper, we aim to extend application of the discontinuous 
Galerkin method to some inverse source problems and for this purpose, we 
present a fully-discrete local discontinuous Galerkin method for solving the 
inverse source problem of the time-fractional diffusion equation. In the pro- 
posed fully-discrete method, we apply a LDG method for the space variable 
and the time-fractional derivative is discretized by using a backward differ- 
ence scheme. The rest of the paper is organized as follows. In Section 2, 
some preliminaries are prepared. Construction of the proposed method for 
the inverse source problem (1)-(3) as well as stability and convergence the- 
orems are dealt with in Section 3. Section 4 is devoted to some numerical 
experiments to illustrate the accuracy and capability of the method. Finally, 
the paper is concluded with a brief conclusion. 


2 Preliminaries 


In this section, some notations are defined and some auxiliary results are 
prepared. At first, we decompose interval [0,1] to some cells (subintervals) 
as follows 

O=71<a3 <--+<ayyi =], 


and set I; = [z;_1, 2542], with the cell lengths Ax; = Ujy1 — @;_1, for 

j =1,...,N and h = maxi<j<n Az;. We denote by ura and ura the 
one 2 2 

values of u at x;,1, from the right cell Jj41 and from the left cell J;. [u)j42 


t+ ul 1, that is the jump of u at the cell interfaces. 
I+3s Its 

For any integer k, we define the piecewise-polynomial space V;* as the space 
of polynomials of degree up to k in each cell J;, thus 


is used to denote wu 


VE = {ve L701]: v|7, €P*(Zj), f=1,...,N}. 


In order to investigate the convergence of the method, we need to define 
projections P and P* as follows 


i (Pu(x) —w(x))v(x)de =0, VYoe P(;), Vi, 


I; 
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i. (P*w(x) — w(x))v(z)dz =0, Vue PP-*(1;), V4, 


j 


Pro(¢t4)= w(x5;51). 


It is well-known that these projections satisfy the following inequality [2,31] 


| 2 |] 4h |] 8 loo the | wo [Inns CARH, 


where w® = Pw —w or w® = P*w — uy, the positive constant C, solely 
depending on w, is independent of h, and 7), denotes the set of boundary 
points of all cells J;. In the following, we use C' to denote a positive constant 
which is independent of h and may have a different value in each occurrence. 
vz denotes the piecewise derivative with respect to x, and the norm || - || 
denotes the usual norm of the space L?{0, 1]. 


3 Construction of the method 


In this section, we construct a numerical scheme for solving problem (1)- 
(3). Let M be a positive integer, At = T/M be the time step size, and 
t, = nAt, n=0,1,...,M denote the time mesh points. An approximation 
to time-fractional derivative (2) can be obtained by a simple quadrature 
formula given as [9,11, 12], 


DP u(2,tn) = ao = p, ilertn—s) <li) O((At)2-%), (5) 


where b; = (¢+1)!~% —i!*. 
Following the LDG regime [6,30], we must rewrite (1) as the following first- 
order system of equations 


q=Us,  Dfu(x,t) — de = f(x)p(t). (6) 


Let u?,gr? € V;* be the approximation of u(-,tn),¢(-,tn) respectively, and 
p” = p(tn),g” = g(tn). Using (5), we establish the necessary weak forms cor- 
responding to (6) and then define a fully-discrete local discontinuous Galerkin 
scheme as follows: find u?,q? € V;*, such that for all test functions v, w € V/*, 
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N 


foiear +3 | f ghesde— She )y44 - Ghe);-) | = 
Q Q jal 
n—-1 ; 
sp" | f(a)vdx + S7(bi-1 - b) | Un ut bn f ubudz, 
Q rar Q Q 

N 

| g;, wdx +f Up, Wedx — So (aw) 544 — (tw );_2) = 0, 
Q Q : 
j=l 


UR (x*) = 9", 

(7) 
where 2 = [0,1] and 6 = (At)°T(2 — a) and without lose of generality we 
assume that x* is a grid point. The “hat” terms in (7) in the cell boundary 
terms from integration by parts are the so-called “numerical fluxes”, which 
are single valued functions defined on the interfaces and should be selected 
carefully for ensuring the numerical stability of the scheme. The choice for 
the numerical fluxes is not unique and among several choices, we can take 
an = (ub) and gt-= (ge)* or ae = (ut )* and @? = (q?)~.. In fact, it is 
important to take i? and gq? from opposite sides [3, 4]. 

Before investigating theoretical aspects of the method, we are going to explain 
details of the method somewhat more. We set 


up (x) = upn(a, nAt) = Yara iF g(x) = an(x, nAt) = Sate 


where N, is the total number of basis functions and 


= ( [tate (e)ae vee [ ieyem,(oree) 


® = (®)(a*) --- ®y (2*)), Z=(0--- 0). 
Setting 6” = (d7 «+» 6h)” and 7” = (77° --- 7m)", scheme (7) leads to the 


following iteration scheme 


n—1 


ky16" + Kioy” = Bp" FP + Sob; _1 — bj) K226"™, 


1=1 


Ko6" pe Koay” = 0, 
05" = g”, 


where n = 1,...,M and 


(Kii)ir = | ©)(x)®,(x)dz, Ko = Ky, 
Q 
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(Kia)ir = af ©;(2)(©,(2))edx — ay (®i(e 5, .)®/ 054) 
— ©(2; 1) ®,r(aF_1) ; 
N 
(Rare = f(a) ()ate-Do (ila) 3) (05, )-le} Jef) 


For solving the direct problem, we have 


n-1 
6” bp’ r+ S (b;-1 — b;)K116"~* Ky, Ky2 
M = , M= , (8 
(S) i=l 0 Ko, Ky (8) 


where matrix Ky; is nonsingular and block diagonal which every block is a 
k xk (k is degree of basis polynomials) matrix. Using the Schur complement, 
linear system (8) has a unique solution iff det(Ki1 — Ki2K{' Ka1) #0. For 
solving the inverse problem, we have 


n—-1 


ky, Ki2 —BF 6” So (bi-1 — b;)Ky16"* 
Ko Ky, ZT yo | =| 1 , (9) 
’ Z 0 p” 0 


g” 
where we put g%’ instead of g” since data are usually obtained by measurement 
tools and have some noises. By solving the nonsymmetric linear system (9) 
using the BiCGStab method, we can obtain uj and p”. 

Just for convenience and without lose of generality, we deal with the case 
g = 0 in the theoretical analysis. In order to examine the stability property 
of the scheme (7), we express following result. 


Theorem 3. Assume that the second derivative of u at x = x* is bounded 
and f is a continuous function on [0,1]. For periodic or compactly supported 
boundary conditions, fully-discrete LDG scheme (7) is unconditionally stable, 
and the numerical solution uj, satisfies 


unl <llugli+s6, n=1,...,M, 


where « is a constant depending on B, f and Uz, at x = x*. 


Proof. We can rewrite scheme (7) as follows 
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N 


[wieder +6 [seeds — SRO Ving - Gey) 
Q 


j=l 
N 
— SS (Gwe )j 44 - (tw) 52) +f gq, wdx +f Up, We dx 
a1 Q Q 
= Bp” | f(x)vdx + So (bi-1 _ bs) f ur—*udx + ina f ul vda, 
: Q Q 
With taking test functions v = u?, w = Gq, for periodic or compactly 
supported boundary conditions, we can obtain 
N 


Bf ated ~ S(aho)j43 - @o*),-4) ) + f ufwade 
Q Q 


j=l 


then 


n-1 
i up, vdx +f gq, wdx = So (bi-1 —b; df un tutda + bn— 7 upunda 
Q Q fy Q 


i=1 
+ sp" | f(x)upda. 
Q 
(10) 
For n = 1 and using Eq. (4), we can get 


Il ub |? +6 ll at [P= 5 ered 26 Jubde 


Q 
rhe te 
<5 ((neen sort wh): 1) Ishi), 
at abas nap 


Using a LDG method for solving an inverse source problem ... 123 


therefore 
I| we WS up I] +6. 


Next, we suppose the following inequalities hold 
ut [<u te, m=1,...,K. 
For n = kK +1, in Eq. (10), and using 


I 


So (bi-1 —b;) + by = 1, 


i=l 


we can obtain 
i 
I] unt? SI] wp I] +8. 


In order to examine the convergence of the scheme (7), we express following 
result. 


Theorem 4. Let u(-,tn) be the exact solution of the problem (1)-(3), which is 
sufficiently smooth with bounded derivatives, and uj, be the numerical solution 
of the fully-discrete LDG scheme (7). There holds the following error estimate 


lal ae SOG (AR? 4 (ADRs 4 e(Ab™), 


where C' is a constant depending ona, u, andT, and c is a constant depending 
on f and Ugg, at x= x*. 


Proof. Obviously for all test functions v,w € V;*, we have 


N 
U(X, tn vax U(X, tn )Wedx — u(z,tn)w);.1 — (u(z,tn)wt);_1 
[vlestayoae + f ula, t)urd Due tabu Yogg = olertad G4) 
N 
+8 | f alestn)ende—S((ale tn) ig — ales ta)o) 9) 
— LO =) f alert ded Onan fale. tadede 
+ [ aertayunde +6 f "(ayuda Splesta) f Fleyvde =o. 


(11) 
Subtracting equation (7) from (11), we can obtain the error equation 
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N 
[ eivae+ J eavede — Y(eRY aay - cerry) 


j=l 
N 
+f jude + ey W,dx — Si(((eR) Ww) 44 - ((e%) w*) 2) 
Q Q jal 
n—-1 
boa | ev udx — So (bi-1 - bs) | en 'uda +6 [ y" (x)udx 
Q = Q Q 


—8p(e.tn) f fa)vde = : 


Cy u(x, tn) ry: Uy = Pe, hs, (P~u(z, tn) St u(z,tn)), (12) 
eG = Qe ty) dr z= Pe, ~~ (Pq(a, tn) _ qa te))i 


Using Eq. (4), we have 


j=1 
N 
+f jude + | en wedx — e”)~w) 41 — ((e") wt), 3 
f egude + f LUM yy EOD ay 
n—-1 
boa | cuode — Y>(bi-1 ~ bi) [ etude +8 [ 9" (a)vdx 
Q mr Q Q 


Q f 


Using Eq. (12), the error equation (13) can be rewritten as follows 


N 
[Pretvae +B (/ Pet unde — )(((Pet)tv-) 44 — (aye) 


j=1 
N 


+ e” wdx er’ w,dxr — er yw) ~e”)~wt), 4a 
[Peiwdes ff Pretunds — (Pret w)j4g ~ (Preey-w),-4) 


j=l 
n-1 
+8 | J salto" ,teude = ee P~ eo vdz + S > (bi-1 - bs) | P~e?'vdz 
af Q oa Q 


+ f (Prue, tn) — u(x, tn))vdx + B ( [atest “qe ta) gar 


J~3 
j=1 


N 
— (Pala, tr) — (a, tn)) TO) jaa — (Pale, tn) — g(a, tn) Fv) 52 } 


+ f (Pale.ta) — q(a,tn))wdx + [@-ueta) — u(x, tn))wrdx 
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= SS (((Po ula, tn) — U(2,tn)) Ww) 54a — (Pula, tn) — wits ty) ) Waa) 


-6 [ 9" (a)uda — tna f (Pula. to) — u(x, to))vdx 


— So (bi-1 — 8) [e-waste — u(x, tp_;))vdz. 


i=l 


With taking the test functions v = u?, w = Gq}, for periodic or compactly 
supported boundary conditions, we derive 


[wren acta f (Pej) ae +a Fittest" sty Pega 
Q 
n-1 
= tna f PoegPmende + (0 i 1-0) [ Pren'Pm ends 
Q 


ign 1 


_B pan (Petar + 00 ( (Pa(z, tn) — a(,tn))* [P-es]) ,_ 3 


2 


+ Lemus ta) —a(@yta))P edt —bn-1 [e-waste — u(x, to))P~ et dx 
~ So (bi-1 b) fe u(x ty i) (x be ;))Po et dx 


For n = 1, we have 


1 @2) ae +8 f ( (Pej) * da ref 4 ~Uge(x*,t1)P e,,dx 


= eo. u(a,to) — u(x, to))P eda 
Q 


3 f Pep ede + BY. ((Palest) ~aee.4))* Pre]), 


j=l is 
- [ Prue,t) —u(z,t,))P~ ede. 
Q 
Using the following facts 
1 
|| Pe |< ChA*, ab < ea? + ee (14) 
€ 


we obtain 


|| Poex (I? +6 ll Peg II’ (| Pen Il +8 I] v(x) || + | Po ule, t1) — u(a, t) | 
+ || P u(x, to) — u(a, to) |] +¢8) || Poe: || 
N 


i 
2 
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< C(h*t! + (At)? + (At) tht? + c(At)®)? + || Poel |I? 
N 


pee 
+BeS> [P eu) j—3 : 
j=l 
If we choose € very small, we conclude that 


|| P- ed ||? +8 || Pel [?< C(n*t! + (At)? + (At) FAPt2 4 c(At)*)?. 
Now, we suppose the following inequalities hold 
|P- er || < C(AP+1 + (At)? + (At) 2 h*+2 + c(At)®%), m=1,2,...,1. 


We need to prove ||P~el*+1|| < C(hk+! + (At)? + (At)th*+2 + c(At)?). 
Letting n =!+1, we have 


[@reyac+e f (Pelt) ae +6 Frttoa( 2" trea )P ey de 
Q Q 
l 


= | P-hPé elt tda + 5" (bi- 1-0) [ prep eltae 
Q 


= ay 


-6 | Yer = Mae +83> (( (Pa(a, ti41) — a(a, tiy1))* ected oa 


j=l 


NIF 


+ [ Poulet) — u(x, tiz1))P ett de 
Q 


by f (Pr u(x, to) —u(a,to))Po elt da 
Q 


i 
= So (bi-1 = bi) [ (Po u(2, tiz1-i) — u(x, ti41—-1))P~ ett de. 
i=1 7 


Then by using (14) and 


we can obtain 
|| Po ett? |? +6 || Pes |? < (| Poe? || + || Po u(a, 41) — u(2, try.) | 
+8 || y'(x) || +b; || Po u(a, to) — u(2, to) || +¢8) || Poet || 
+1 = 6) | Poe Poet | Sr dey? 
= 


iS N 

B 
ae ae (Pa(2, ti41) — g(a, tit1))*) 54 
j= 
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< Cyb (n't! + (At)? + (At) Fh*+2 + c(At)®) || Poel? | 

> [Pe], a + (bi Sb Cathe ths (AD (AD 2 AP 
+o(At)®) (Peis Cc + (At)? + (At) #nk+2 + c(At)*)? 

be || Po elt! |]? 3 [Po elt]? 


j-3¢ 


d=! 


Choosing a small ¢, we derive 


|P~elt2|| < C(hR+1 + (At)? + (At) Akt? 4+ c(At)®). 


4 Numerical examples 


In this section, we carry out some numerical tests to confirm theoretical re- 
sults and to investigate the efficiency of the proposed method. The maximum 
time is J’ = 1 otherwise it will be specified. The space and time step sizes 
are h = 1/N and At = T/M, respectively. We use the relative root mean 
square error, i.e., 


M M 1/2 
e(p) = (Sous - and? Zoot) (15) 


n=1 n=1 


for checking the accuracy of the numerical solutions. In all of tests, we take 
k = 2, i.e., we consider piecewise polynomials of degree two as the basis 
functions in the LDG regime. For dealing with the sensitivity of the solution 
with respect to the data, we use the following noisy data 


g°(tn) = g(tn)(1+érnd(n)), n=0,1,..., 


where g is the exact data and rnd(n) is a random number uniformly dis- 
tributed in [—1, 1] and the magnitude 6 indicates a relative noise level. 


Example 1 We consider the inverse source problem of the time-fractional 
diffusion equation (1)-(3) with the exact solution u(z,t) = e~' cos(27zx). Set- 
ting At very small and using the usual L? and L® error norms, we show 
in Table 1 that the order of convergence of the proposed method is about 
three as we expected (according to the obtained error estimate since k = 2). 
Since the exact p of this problem is not accessible, in Fig 1. we show pj} for 
a = 0.5,2* = 0.75 and N = 8,16. The errors in L?-norm and L®-norm for 
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piecewise P*,k = 1,2,3 polynomials for a = 0.1,2* = 0.25 are presented in 
Fig 2. 


Ne8, o=0.5, x'=0.75 Ne16, a=0.5, x'=0.75 


45 45 

40 40 

35 35 

30 30 

25 25 

= 2 

20 20 

15 15 

10 10 

5 5 

% 0.2 0.4 ; 0.6 0.8 1 % 0.2 0.4 : 0.6 0.8 1 
Figure 1: Numerical approximations to p for Example 1 for N = 8 (left) and N = 16 


(right). 


-2 -2 
—4 4 4 —4Ff J 
-6+ 4 -6t 4 
S -s} { BS —gt 4 
© rr) 
a ov 
B -10} { B-10}+ | 
12+ 4 12} | 
—-14b es ks J —14} 4 
—@— k=2 
—<— k=3 
16 16 
-4 = 2 1 4 3 = =i 
log(h) log(h) 


Figure 2: log(error) as a function of log(h) for a = 0.1,2* = 0.25 when using piecewise 
P*,k =1,2,3 polynomials for Example 1. 


Example 2. Let the exact solution for problem (1)-(3) be u(a,t) = 
#? sin (Za). Therefore, (x) = u(x,0) = 0, ko(t) = u(0,t) = 0, k(t) = 


u(1, t) = ae f(x) — sin ($2), and p(t) — 
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Table 1: Accuracy test for Example 1 for different a with 2* = 0.5. 


N | L? error | Order | L™ error | Order 
a=0.3 | 5 | 0.002067 - 0.002780 - 
10 | 0.000280 2.9 0.000376 2.9 
15 | 0.000079 3.1 0.000112 3.0 
20 | 0.000034 2.9 0.000047 3.0 
a=0.5 | 5 | 0.002068 - 0.002781 - 

10 | 0.000280 2.9 0.000376 2.9 
15 | 0.000079 3.1 0.000112 3.0 
20 | 0.000034 2.9 0.000047 3.0 
a=0.7 | 5 | 0.002069 - 0.002783 - 

10 | 0.000280 2.9 0.000377 2.9 
15 | 0.000079 3.1 0.000112 3.0 
20 | 0.000034 2.9 0.000047 3.0 
a=1 5 | 0.002070 - 0.002784 - 

10 | 0.000281 2.9 0.000377 2.9 
15 | 0.000079 3.1 0.000112 3.0 
20 | 0.000034 2.9 0.000047 3.0 


x* = 0.5, then g(t) = sin(7)t?. p? for a = 0.5 and a = 0.95 with various noise 
levels 6 = 5%, 10%, 15% are plotted in Fig. 3. Corresponding relative root 
mean square errors are €(p) = 8.5762 x 10~°, 8.6406 x 10~°, 8.7104 x 10~° 
for a = 0.5 and e(p) = 2.030156 x 103, 3.077799 x 10-3, 4.153917 x 10~8 for 
a = 0.95. In Table 2, we compare the relative root mean square errors for the 


X*=0.5, a=0.5, N=50, M=200 X*=0.5, a=0.95, N=50, M=200 


Figure 3: Numerical approximations to p for Example 2. 
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proposed method in [28] (€1(p) for no regularization method and ¢€2(p) with 
a regularization method) with the proposed LDG method (¢3(p)). In Table 3 
the relative errors for different x* with N = 50 are compared. The proposed 
method could generate more satisfactory results without any regularization 
method. To verify the role of the final time T, we depict pj for a = 0.95 
with T = 100, 10000 and 6 = 5%, 10%, 15% in Fig. 4. We can see that the 
dependency of the results to final time T is almost unimportant, even when 
T is very large, i.e., T = 10000. Without any regularization method, the 
trace of the loss of stability does not appear. 


Table 2: Comparison approximation solutions of Example 2 with the literature for various 


a. 

a=0.1 a=0.3 a=0.5 a=0.7 a=0.9 
5=5% | e1(p) 0.0773 0.0880 0.1230 0.2651 0.8079 
€2(p) 0.033 0.0362 0.0393 0.0413 0.0431 

e3(p) | 2.162x10~-° | 1.9319x 107° | 8.5698 x 107° | 3.33882 x 107* | 0.001419157 
5=10% | e1(p) 0.1545 0.1738 0.2436 0.5296 1.6159 
€2(p) 0.0634 0.0633 0.0662 0.0721 0.0812 

€3(p) | 2.1600 x 107° | 1.9334 x 1075 | 8.6232 x 107° | 3.51023 x 1074 | 0.002105562 
5=15% | €1(p) 0.2322 0.2611 0.3658 0.7948 2.4241 
€2(p) 0.0952 0.0943 0.0948 0.0721 0.1248 

e3(p) | 2.1660 x 10~° | 1.9356 x 107° | 8.6557 x 10~° | 3.96662 x 10~* | 0.002444287 


Table 3: Comparison approximation solutions of Example 2 with the literature for various 


c 


x* | e1(p) | é2(p) €3(p) 
0.1 | 1.0010 | 0.1801 | 1.024501 x 107? 
0.2 | 0.9047 | 0.0661 | 8.97953 x 10~* 
0.3 | 0.9197 | 0.0767 | 8.03529 x 10~° 
0.4 | 0.9323 | 0.0694 | 7.96885 x 10~* 
0.5 | 0.9222 | 0.0839 | 7.50627 x 10~* 
0.6 | 0.9071 | 0.0761 | 7.92723 x 10~* 
0.7 | 0.9147 | 0.1026 | 7.78528 x 10~° 
0.8 | 0.9256 | 0.0974 | 7.55676 x 10~° 


(1)-(3), with 
0, f(x) = 27, 


Example 3. We test a none-smooth problem corresponding to 
(x) = u(x, 0) = sin(2rax), ko(t) = u(0,t) = 0, k(t) = u(t) = 
and 
@={"re t € [0, 0.5), 
—2t+2+a,t € (0.5, 1]. 
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Figure 4: Numerical approximations to p for Example 2 for T = 100 (left) and T = 10000 
(right). 


Since the exact solution of this problem is not accessible, we first solve a 
direct problem using a suitable LDG method to obtain the input data g 
then we solve the inverse problem using our method. In [28], the direct 
problem has been solved by using an implicit finite difference (FD) method 
but since p is not a smooth function we expect to have a none-smooth solution 
and we decided to solve the direct problem using a LDG method too. We 
have to point out that since in our method we face the sparse systems, the 
computational complexity of both methods, i.e., FD and LDG is almost equal. 
py for a = 0.5, 0.95 with noise levels 6 = 5%, 10%, 15% are presented in 
Fig 5. Without applying any regularization methods, our results are in good 
agreement with the results of [28]. The corresponding relative root mean 
square errors are €(p) = 3.8300 x 1077, 6.2300 x 1077, 9.4800 x 107” for 
a = 0.5 and e(p) = 1.03320 x 1074, 2.33279 x 1074, 1.03320 x 1074 for 
a = 0.95, which are considerably better than reported in [28]. 


Example 4. We test a discontinuous problem corresponding to (1)-(3), with 
(x) = u(x, 0) = sin(2ma), ko(t) = u(0, t) = 0, kx (t) = u(1,t) = 0, f(a) = 2”, 
and 
(t) = { 1, t € [0.25, 0.75], 
~ | 0, t € [0, 0.25) U (0.75, 1]. 


Since the exact solution of this problem is not accessible, we first solve a 
direct problem using a suitable LDG method to obtain the input data g 
then we solve the inverse problem using our method. In [28], the direct 
problem has been solved by using an implicit finite difference method but 
since p is not a continuous function we expect to have a discontinuous solution 
and we decided to solve the direct problem using a LDG method too. We 
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Figure 5: Numerical approximations to p for Example 3 with various a. 


have to point out that since in our method we face the sparse systems, the 
computational complexity of both methods, i.e., FD and LDG is almost equal. 
pp for a = 0.5, 0.95 with noise levels 6 = 5%, 10%, 15% are plotted in Fig 
6. Without applying any regularization methods, our results are in good 
agreement with the results of [28]. The corresponding relative root mean 
square errors are e(p) = 4.1200 x 1077, 6.2300 x 1077, 8.2700 x 107” for 
a = 0.5 and e(p) = 8.8388 x 1075, 1.61535 x 1074, 2.58527 x 107+ for 
a = 0.95, which are considerably better than reported in [28]. 


X*=0.5, «=0.5, N=200, M=200 X*=0.5, «=0.95, N=200, M=200 
1.2 1.2 
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Figure 6: Numerical approximations to p for Example 4 with various a. 


Using a LDG method for solving an inverse source problem ... 133 


5 Conclusion 


In this paper, an inverse source problem for the time-fractional diffusion equa- 
tion was solved numerically by using a local discontinuous Galerkin method. 
In fact, we could extend a fully-discrete LDG finite element method for solv- 
ing a class of time-fractional inverse problem. By applying this method with- 
out using any regularization methods, we could obtain stable and accurate 
numerical approximations to the time-dependent source term using an addi- 
tional data in an interior measurement location. The numerical stability and 
convergence of the proposed method have been investigated and theoretically 
proven. Various numerical examples with smooth or none-smooth data and 
maybe solutions have been verified to demonstrate the effectiveness and ro- 
bustness of the proposed method. This outstanding and promising method 
can be further applied to another one-dimensional or higher dimensional in- 
verse problems which can be considered for the future works. 
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